---
title: "IH Models (Paper 2)"
author: "Samuel Johnston"
date: "07/11/2020"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

## R Markdown

This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see <http://rmarkdown.rstudio.com>.

When you click the **Knit** button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:

```{r cars}
summary(cars)
```

## Including Plots

You can also embed plots, for example:

```{r pressure, echo=FALSE}
plot(pressure)
```

Note that the `echo = FALSE` parameter was added to the code chunk to prevent printing of the R code that generated the plot.

```{r}
#Packages used
library(dplyr)
library(plm)
library(stargazer)
library(multiwayvcov)
library(lmtest)
library(sandwich)
library(Rmisc)
library(ggplot2)
library(ggeffects)
```


```{r}
#LOADING DATA AND RECODING VARIABLES
library(foreign)
ihdata <- read.csv("/Users/samueljohnston/Documents/R Code/Thesis/Internal Homogenisation Dataset.csv")
descriptive <- read.csv("/Users/samueljohnston/Documents/R Code/Thesis/Descriptive Data.csv")
allotherparties <- read.csv("/Users/samueljohnston/Documents/R Code/Thesis/All Other Parties IH.csv")

#In order to recode the variables, in terms of simplifying variable names and ensuring they are numeric:
ihdata$Centre_left_decent_sal <- as.numeric(ihdata$Centre_left_decent_sal)
ihdata$GDP_per_capita <- as.numeric(ihdata$GDP_per_capita)

allotherparties$Centre.left.decent..sal. <- as.numeric(allotherparties$Centre.left.decent..sal.)
allotherparties$GDP.per.capita <- as.numeric(allotherparties$GDP.per.capita)


#To transform the IH variables from a proportion to a percentage:
ihdata$IH_Position <- ihdata$IH_Position*100
ihdata$Pro_IH <- ihdata$Pro_IH*100
descriptive$IH.Position <- descriptive$IH.Position*100
allotherparties$IH.Position <- allotherparties$IH.Position*100

#Creating a descriptive plot (Figure 1):
plot <- ggplot(descriptive, aes(x = IH.Position, y = X5.Years, linetype = factor(EU.Membership))) +
  labs(x = "Relative Saliency of Ethnoregionalism", y = "Five Year Interval", title = "Figure 1: Trends in the Relative Saliency of Ethnoregionalism Over Time,", subtitle = "Contingent on European Union Membership") + 
  scale_linetype_discrete(name = "EU Membership") +
  geom_line() +
  theme_classic()
plot + geom_jitter()
ggsave(file = '/Users/samueljohnston/Documents/R Code/Thesis/Figure 1.eps')
```


```{r}
#To run the models for IH and EU membership (Hypothesis 1), first without the interaction, then with (robust SEs follow each model): 
model1 <- plm(IH_Position ~ EU_Member_Lagged + RAI + Lagged_IH + Lagged_RR_Presence + Lagged_ER_Presence + Established_Party_Distance + Centre_left_decent_sal + Centre_right_decent_sal + International_migrant_stock + ELF + Economic_growth + GDP_per_capita + Unemployment + Disproportionality + parfam + Eastern_European, data = ihdata, model = "within", index = "Election_Year")

#This tests whether the model suffers from heteroskedasticity (Breusch-Pagan test), and it shows that we do, so we need to use robust standard errors.
bptest(IH_Position ~ EU_Member_Lagged + RAI + Lagged_IH + Lagged_RR_Presence + Lagged_ER_Presence + Established_Party_Distance + Centre_left_decent_sal + Centre_right_decent_sal + International_migrant_stock + ELF + Economic_growth + GDP_per_capita + Unemployment + Disproportionality + parfam + Eastern_European + factor(Election.Year), data = ihdata, studentize = F) 

cov.fit1 <- vcovHC(model1, type = "HC0")
rob.std.err1 <- sqrt(diag(cov.fit1))

model2 <- plm(IH_Position ~ EU_Member_Lagged*RAI + Lagged_IH + Lagged_RR_Presence + Lagged_ER_Presence + Established_Party_Distance + Centre_left_decent_sal + Centre_right_decent_sal + International_migrant_stock + ELF + Economic_growth + GDP_per_capita + Unemployment + Disproportionality + parfam + Eastern_European, data = ihdata, model = "within", index = "Election_Year")

cov.fit2 <- vcovHC(model2, type = "HC0")
rob.std.err2 <- sqrt(diag(cov.fit2))


stargazer(model1, model2, se = list(rob.std.err1, rob.std.err2), dep.var.labels = c("Relative Saliency of Ethnoregionalism"), covariate.labels = c("EU Membership", "Regionalisation", "Lagged Dependent Variable", "Lagged Radical Right Presence", "Lagged Ethnoregional Presence", "Established Party Convergence", "Centre-Left Decentralisation Saliency", "Centre-Right Decentralisation Saliency", "Immigration", "Ethnolinguistic Fractionalisation", "GDP Growth", "GDP per capita", "Unemployment", "Disproportionality", "Party Family", "Eastern European", "EU*Regionalisation Interaction"),  title = "Table 1: Determinants of the Relative Saliency of Ethnoregionalism for 33 European Countries, 1980-2019", no.space = TRUE, style = "apsr", star.char = c("*", "**"), star.cutoffs = c(.05, .01), notes = c("Both models use OLS and robust standard errors", "*p < .05; **p < .01"), notes.append = FALSE, type = "text", out = "tableIH.txt")

summary(ihdata$RAI, na.rm = TRUE)
sd(ihdata$RAI, na.rm = TRUE)
```

```{r}
#Does EU Membership have an effect if we just look at cultural protectionism?  
model1B <- plm(CP_Saliency ~ EU_Member_Lagged + RAI + Lagged_CP_Saliency + Lagged_RR_Presence + Established_Party_Distance + Centre_left_cult_prot_sal + Centre_right_cult_prot_sal + International_migrant_stock + ELF + Economic_growth + GDP_per_capita + Unemployment + Disproportionality + parfam + Eastern_European, data = ihdata, model = "within", index = "Election_Year")

cov.fit1B <- vcovHC(model1B, type = "HC0")
rob.std.err1B <- sqrt(diag(cov.fit1B))

model1C <- plm(CP_Saliency ~ EU_Member_Lagged*RAI + Lagged_CP_Saliency + Lagged_RR_Presence + Established_Party_Distance + Centre_left_cult_prot_sal + Centre_right_cult_prot_sal + International_migrant_stock + ELF + Economic_growth + GDP_per_capita + Unemployment + Disproportionality + parfam + Eastern_European, data = ihdata, model = "within", index = "Election_Year")

cov.fit1C <- vcovHC(model1C, type = "HC0")
rob.std.err1C <- sqrt(diag(cov.fit1C))


stargazer(model1B, model1C, se = list(rob.std.err1B, rob.std.err1C), dep.var.labels = c("Cultural Protectionism Saliency"), covariate.labels = c("EU Membership", "Regionalisation", "Lagged Cultural Protectionism Saliency", "Lagged Radical Right Presence", "Established Party Convergence", "Centre-Left Cultural Protectionism Saliency", "Centre-Right Cultural Protectionism Saliency", "Immigration", "Ethnolinguistic Fractionalisation", "GDP Growth", "GDP per capita", "Unemployment", "Disproportionality", "Party Family", "Eastern European", "EU Membership*Regionalisation Interaction"), title = "Table 2: Determinants of the Saliency of Cultural Protectionism for 33 European Countries, 1980-2019", no.space = TRUE, style = "apsr", star.char = c("*", "**"), star.cutoffs = c(.05, .01), notes = c("Both models use OLS and robust standard errors", "*p < .05; **p < .01"), notes.append = FALSE, type = "text", out = "table2.txt")
```

```{r}
#APPENDIX.
#DESCRIPTIVE PLOT (Figure 3)
ihdataB <- ihdata %>% filter(ihdata$Country %in% c("Belgium", "Finland", "Italy", "Norway", "Poland", "Romania", "Slovakia", "Spain", "Switzerland"))
ihdataC <- ihdataB %>% filter(ihdataB$Party.Abbrev %in% c("VB", "RKP/SFP", "LN", "FrP", "PiS", "UDMR/RMDSz", "SNS", "SMK/MKP", "PNV/EAJ", "ERC", "SVP/UDC"))

plot <- ggplot(ihdataC, aes(x = IH_Position, y = Election_Year, group = factor(Party))) +
  labs(x = "Relative Saliency of Ethnoregionalism", y = "Year", title = "Figure 3: Trends in the Relative Saliency of Ethnoregionalism Over Time,", subtitle = "for a Selection of Parties") + 
  geom_point(aes(colour = Party))
ggsave(file = '/Users/samueljohnston/Documents/R Code/Thesis/Figure 3.eps')


#ROBUSTNESS CHECK 1: To check if the results are driven merely by ethnoregional parties, I subset the data to replicate the previous regressions only for radical right parties: 
ihdata3 <- ihdata %>% filter(! ihdata$parfam == 0)

model1D <- plm(IH_Position ~ EU_Member_Lagged + RAI + Lagged_IH + Lagged_RR_Presence + Lagged_ER_Presence + Established_Party_Distance + Centre_left_decent_sal + Centre_right_decent_sal + International_migrant_stock + ELF + Economic_growth + GDP_per_capita + Unemployment + Disproportionality + Eastern_European, data = ihdata3, model = "within", index = "Election_Year")

cov.fit1D <- vcovHC(model1D, type = "HC0")
rob.std.err1D <- sqrt(diag(cov.fit1D))

model1E <- plm(IH_Position ~ EU_Member_Lagged*RAI + Lagged_IH + Lagged_RR_Presence + Lagged_ER_Presence + Established_Party_Distance + Centre_left_decent_sal + Centre_right_decent_sal + International_migrant_stock + ELF + Economic_growth + GDP_per_capita + Unemployment + Disproportionality + Eastern_European, data = ihdata3, model = "within", index = "Election_Year")

cov.fit1E <- vcovHC(model1E, type = "HC0")
rob.std.err1E <- sqrt(diag(cov.fit1E))

stargazer(model1D, model1E, se = list(rob.std.err1D, rob.std.err1E), dep.var.labels = c("Internal Homogenisation Position"), covariate.labels = c("EU Membership", "Regionalisation", "Lagged Dependent Variable", "Lagged Radical Right Presence", "Lagged Ethnoregional Presence", "Established Party Convergence", "Centre-Left Decentralisation Saliency", "Centre-Right Decentralisation Saliency", "Immigration", "Ethnolinguistic Fractionalisation", "GDP Growth", "GDP per capita", "Unemployment", "Disproportionality", "Eastern European", "EU Membership*Regionalisation Interaction"), title = "Table 3: Determinants of Internal Homogenisation Position for RRPs in 33 European Countries, 1980-2019", no.space = TRUE, style = "apsr", star.char = c("*", "**"), star.cutoffs = c(.05, .01), notes = c("Both models use OLS and robust standard errors", "*p < .05; **p < .01"), notes.append = FALSE, type = "text", out = "table3.txt")
```

```{r}
#ROBUSTNESS CHECK 2: REPLACING EASTERN EUROPE WITH OTHER CONTROLS (Austro-Hungarian Empire, Ottoman Empire, Yugoslavia)
model3 <- plm(IH_Position ~ EU_Member_Lagged*RAI + Lagged_IH + Lagged_RR_Presence + Lagged_ER_Presence + Established_Party_Distance + Centre_left_decent_sal + Centre_right_decent_sal + International_migrant_stock + ELF + Economic_growth + GDP_per_capita + Unemployment + Disproportionality + parfam + Austro_Hungarian_Empire, data = ihdata, model = "within", index = "Election_Year")

cov.fit3 <- vcovHC(model3, type = "HC0")
rob.std.err3 <- sqrt(diag(cov.fit3))

model4 <- plm(IH_Position ~ EU_Member_Lagged*RAI + Lagged_IH + Lagged_RR_Presence + Lagged_ER_Presence + Established_Party_Distance + Centre_left_decent_sal + Centre_right_decent_sal + International_migrant_stock + ELF + Economic_growth + GDP_per_capita + Unemployment + Disproportionality + parfam + Ottoman_Empire, data = ihdata, model = "within", index = "Election_Year")

cov.fit4 <- vcovHC(model4, type = "HC0")
rob.std.err4 <- sqrt(diag(cov.fit4))

model5 <- plm(IH_Position ~ EU_Member_Lagged*RAI + Lagged_IH + Lagged_RR_Presence + Lagged_ER_Presence + Established_Party_Distance + Centre_left_decent_sal + Centre_right_decent_sal + International_migrant_stock + ELF + Economic_growth + GDP_per_capita + Unemployment + Disproportionality + parfam + Yugoslavia, data = ihdata, model = "within", index = "Election_Year")

cov.fit5 <- vcovHC(model5, type = "HC0")
rob.std.err5 <- sqrt(diag(cov.fit5))

stargazer(model3, model4, model5, se = list(rob.std.err3, rob.std.err4, rob.std.err5), dep.var.labels = c("Relative Saliency of Ethnoregionalism"), covariate.labels = c("EU Membership", "Regionalisation", "Lagged Dependent Variable", "Lagged Radical Right Presence", "Lagged Ethnoregional Presence", "Established Party Convergence", "Centre-Left Decentralisation Saliency", "Centre-Right Decentralisation Saliency", "Immigration", "Ethnolinguistic Fractionalisation", "GDP Growth", "GDP per capita", "Unemployment", "Disproportionality", "Party Family", "Austro-Hungarian Empire", "Ottoman Empire", "Yugoslavia", "EU*Regionalisation Interaction"), title = "Table 4: Determinants of the Relative Saliency of Ethnoregionalism for Thirty-Three European Countries, 1980-2019", no.space = TRUE, style = "apsr", star.char = c("*", "**"), star.cutoffs = c(.05, .01), notes = c("All models use OLS and robust standard errors", "*p < .05; **p < .01"), notes.append = FALSE, type = "text", out = "table4.txt")
```


```{r}
#ROBUSTNESS CHECK 3: ALTERNATIVE CONTROLS
model6 <- plm(IH_Position ~ EU_Member_Lagged*RAI + Lagged_IH + Established_Party_Distance + International_migrant_stock + ELF + GDP_per_capita + Unemployment + Disproportionality + parfam + Eastern_European, data = ihdata, model = "within", index = "Election_Year")

cov.fit6 <- vcovHC(model6, type = "HC0")
rob.std.err6 <- sqrt(diag(cov.fit6))

model7 <- plm(IH_Position ~ EU_Member_Lagged*RAI + Lagged_IH + Established_Party_Distance + Centre_left_decent_sal + Centre_right_decent_sal + International_migrant_stock + ELF + GDP_per_capita + Unemployment + Disproportionality + parfam + Eastern_European, data = ihdata, model = "within", index = "Election_Year")

cov.fit7 <- vcovHC(model7, type = "HC0")
rob.std.err7 <- sqrt(diag(cov.fit7))

model8 <- plm(IH_Position ~ EU_Member_Lagged*RAI + Lagged_IH + Established_Party_Distance + Centre_left_decent_sal + Centre_right_decent_sal + International_migrant_stock + ELF + Economic_growth + GDP_per_capita + Unemployment + Disproportionality + parfam + Eastern_European, data = ihdata, model = "within", index = "Election_Year")

cov.fit8 <- vcovHC(model8, type = "HC0")
rob.std.err8 <- sqrt(diag(cov.fit8))


model9 <- plm(IH_Position ~ EU_Member_Lagged*RAI + Lagged_IH + Established_Party_Distance + + Centre_left_decent_sal + Centre_right_decent_sal + International_migrant_stock + ELF + Economic_growth + GDP_per_capita + Unemployment + Disproportionality + parfam2 + Eastern_European, data = ihdata, model = "within", index = "Election_Year")

cov.fit9 <- vcovHC(model9, type = "HC0")
rob.std.err9 <- sqrt(diag(cov.fit9))


stargazer(model6, model7, model8, model9, se = list(rob.std.err6, rob.std.err7, rob.std.err8, rob.std.err9), dep.var.labels = c("Relative Saliency of Ethnoregionalism"), covariate.labels = c("EU Membership", "Regionalisation", "Lagged Dependent Variable", "Established Party Convergence", "Centre-Left Decentralisation Saliency", "Centre-Right Decentralisation Saliency", "Immigration", "Ethnolinguistic Fractionalisation", "GDP Growth", "GDP per capita", "Unemployment", "Disproportionality", "Party Family", "Alternative Party Family", "Eastern European", "EU Membership*Regionalisation Interaction"), title = "Table 5: Determinants of the Relative Saliency of Ethnoregionalism for Thirty-Three European Countries, 1980-2019", no.space = TRUE, star.char = c("*", "**"), star.cutoffs = c(.05, .01), notes = c("All models use OLS and robust standard errors", "*p < .05; **p < .01"), notes.append = FALSE, style = "apsr", type = "text", out = "table5.txt")
```

```{r}
#ROBUSTNESS CHECK 4: EXCLUDING ANTI-CATEGORIES FROM DV AND RATIO DV
#Firstly, we want to ensure that including the anti-nationalist categories from MARPOR (i.e., per302 is pro-centralisation, per602 is anti-national way of life, and per607 is pro-multiculturalism) in the construction of the DV isn't biasing our results. Consequently, I replicate the models where only per301 (pro-decentralisation), per601 (pro the national way of life), and per608 (anti-multiculturalism) are used to construct the DV. The party competition controls are not changed as other parties may help to politicise the issue by taking the opposite position on it, so a wider saliency measure is still informative: 
model10 <- plm(Pro_IH ~ EU_Member_Lagged*RAI + Lagged_Pro_IH + Lagged_RR_Presence + Lagged_ER_Presence + Established_Party_Distance + Centre_left_decent_sal + Centre_right_decent_sal + International_migrant_stock + ELF + Economic_growth + GDP_per_capita + Disproportionality + parfam + Eastern_European, data = ihdata, model = "within", index = "Election_Year")

cov.fit10 <- vcovHC(model10, type = "HC0")
rob.std.err10 <- sqrt(diag(cov.fit10))

#Secondly, in order to ensure that we get the same results if we use a ratio, rather than a proportion: 
model11 <- plm(IH_Ratio ~ EU_Member_Lagged*RAI + Lagged_IH_Ratio + Lagged_RR_Presence + Lagged_ER_Presence + Established_Party_Distance + Centre_left_decent_sal + Centre_right_decent_sal + International_migrant_stock + ELF + Economic_growth + GDP_per_capita + Disproportionality + parfam + Eastern_European, data = ihdata, model = "within", index = "Election_Year")

cov.fit11 <- vcovHC(model11, type = "HC0")
rob.std.err11 <- sqrt(diag(cov.fit11))

stargazer(model10, model11, se = list(rob.std.err10, rob.std.err11), dep.var.labels = c("Pro-Nationalism Categories", "IH Ratio"),  title = "Table 6: Determinants of the Relative Saliency of Ethnoregionalism for Thirty-Three European Countries, 1980-2019", covariate.labels = c("EU Membership", "Regionalisation", "Lagged Pro-Categories", "Lagged Ratio", "Lagged RRP Presence", "Lagged Ethnoregional Presence", "Established Party Convergence", "Centre-Left Decentralisation Saliency", "Centre-Right Decentralisation Saliency", "Immigration", "Ethnolinguistic Fractionalisation", "Economic Growth", "GDP per capita", "Disproportionality", "Party Family", "Eastern European", "EU Membership*Regionalisation Interaction"), no.space = TRUE, star.char = c("*", "**"), star.cutoffs = c(.05, .01), notes = c("All models use OLS and robust standard errors", "*p < .05; **p < .01"), notes.append = FALSE, style = "apsr", type = "text", out = "table6.txt")
```


```{r}
#ROBUSTNESS CHECK 5: ALL OTHER PARTIES
model12 <- plm(IH.Position ~ EU.Member_Lagged*RAI + as.numeric(Lagged.IH.Position) + Lagged.RR.Presence + Lagged.ER.Presence + Established.Convergence + Centre.left.decent..sal. + Centre.right.decent..sal. + International.migrant.stock + ELF + Economic.Growth + Unemployment + Disproportionality + Eastern.European, data = allotherparties, model = "within", index = "Year")

cov.fit12 <- vcovHC(model12, type = "HC0")
rob.std.err12 <- sqrt(diag(cov.fit12))

stargazer(model12, se = list(rob.std.err12), dep.var.labels = c("Relative Saliency of Ethnoregionalism"), covariate.labels = c("EU Membership", "Regionalisation", "Lagged Dependent Variable", "Radical Right Presence", "Ethnoregional Presence", "Established Party Convergence", "Centre-left Decentralisation Saliency", "Centre-right Decentralisation Saliency", "Immigration", "Ethnolinguistic Fractionalisation", "GDP Growth", "Unemployment", "Disproportionality", "Eastern European", "EU Membership*Regionalisation Interaction"), title = "Table 7: Determinants of the Relative Saliency of Ethnoregionalism for Thirty-Three European Countries, 1980-2019", no.space = TRUE, star.char = c("*", "**"), star.cutoffs = c(.05, .01), notes = c("All models use OLS and robust standard errors", "*p < .05; **p < .01"), notes.append = FALSE, style = "apsr", type = "text", out = "table7.txt")
```

